home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-01 / ddj9304.zip / 1993-APR.ZIP / STEREO.ASC < prev    next >
Text File  |  1993-02-17  |  15KB  |  388 lines
  1. _ALGORITHMS FOR STEREOSCOPIC IMAGES_
  2. by Victor Duvanenko and W.E. Robbins
  3.  
  4. [LISTING ONE]
  5.  
  6. display_3D_data( proj )
  7. int proj;
  8. {
  9.    double  Tw[4][4], S[4][4], Td[4][4], Ry[4][4], Per[4][4], tmp[4][4];
  10.    double  Left[4][4], Right[4][4];
  11.    double  e_v,      /* interocular distance mapped back to the view  coord. */
  12.            e_w;      /* interocular distance mapped back to the world coord. */
  13.    printf( "Computing normals for shading.  " );
  14.    compute_normals();      /* and the dot products for each triangle */
  15.    printf( "Done.\n" );
  16.    /* Perspective projection must use three steps: 
  17.       1)  Compute normals and project,
  18.       2)  Divide by W (homogenize)
  19.       3)  Transform into the device coordinates. */
  20.    if ( proj == PERSPECTIVE )
  21.    {
  22.       /* map physical interocular distance into the view port coordinates. */
  23.       e_v = INTEROCULAR_DISTANCE / PIXEL_WIDTH;      /* e_v == pixels */
  24.       /* map from the viewport coordinate system to the world. */
  25.       e_w = e_v / (( x_right_d - x_left_d ) / ( x_right - x_left ));
  26.       set_to_identity( Left,  4 );
  27.       set_to_identity( Right, 4 );
  28.  
  29.       /* Use the Translate, Project, Translate back model. */
  30.  
  31.       /* Create the Left eye transformation matrix. */
  32.       set_to_identity( Tw, 4 );    /* translate the left eye to the origin */
  33.       Tw[3][0] = -e_w / 2.0;
  34.       matrix_mult( Left, 4, 4, Tw, 4, 4, tmp );
  35.       /* Create the perspective projection matrix. */
  36.       set_to_identity( Per, 4 );
  37.       Per[2][3] = 1.0 / proj_plane;          /* 1/d */
  38.       Per[3][3] = 0.0;
  39.       matrix_mult( tmp, 4, 4, Per, 4, 4, Left );
  40.  
  41.       Tw[3][0] = e_w / 2.0;             /* translate back */
  42.       matrix_mult( Left, 4, 4, Tw, 4, 4, tmp );
  43.       copy_matrix( tmp, Left, 4, 4 );
  44.  
  45.       /* Create the Right eye transformation matrix. */
  46.       set_to_identity( Tw, 4 );    /* translate the right eye to the origin */
  47.       Tw[3][0] = e_w / 2.0;
  48.       matrix_mult( Right, 4, 4, Tw, 4, 4, tmp );
  49.       /* Create the perspective projection matrix. */
  50.       set_to_identity( Per, 4 );
  51.       Per[2][3] = 1.0 / proj_plane;          /* 1/d */
  52.       Per[3][3] = 0.0;
  53.       matrix_mult( tmp, 4, 4, Per, 4, 4, Right );
  54.  
  55.       Tw[3][0] = -e_w / 2.0;             /* translate back */
  56.       matrix_mult( Right, 4, 4, Tw, 4, 4, tmp );
  57.       copy_matrix( tmp, Right, 4, 4 );
  58. #if 0
  59.       printf( "Transforming, projecting and homogenizing the image.  " );
  60.       transform_and_homogenize_image( image, tr_image, Tm );
  61.       printf( "Done.\n" );
  62. #endif
  63.    }
  64.    /* Create the world to device transformation matrix. */
  65.    /* Create the translation matrix. Translate only in X and Y. */
  66.    set_to_identity( Tw, 4 );
  67.    Tw[3][0] = -x_left;   Tw[3][1] = -y_bottom;   Tw[3][2] = 0.0;
  68.    /* Create a uniform scale matrix. */
  69.    set_to_identity( S, 4 );
  70.    S[0][0] = ( x_right_d - x_left_d   ) / ( x_right - x_left   );
  71.    S[1][1] = ( y_top_d   - y_bottom_d ) / ( y_top   - y_bottom );
  72.    S[2][2] = ( z_back_d  - z_front_d  ) / ( z_back  - z_front  );
  73.  
  74.    matrix_mult( Tw, 4, 4, S, 4, 4, tmp );
  75.    copy_matrix( tmp, Tw, 4, 4 );
  76.  
  77.    /* Create the translation matrix. Translate only in X and Y. */
  78.    set_to_identity( Td, 4 );
  79.    Td[3][0] = x_left_d;   Td[3][1] = y_bottom_d;   Td[3][2] = 0.0;
  80.    matrix_mult( Tw, 4, 4, Td, 4, 4, tmp );
  81.    copy_matrix( tmp, Tw, 4, 4 );
  82.  
  83.    /* Since device/screen origin on LEX/90 is in upper left, we need to reflect
  84.       Y and translate by screen height to place device origin in bottom left.*/
  85.    set_to_identity( Ry, 4 );
  86.    Ry[1][1] = -1.0;
  87.    matrix_mult( Tw, 4, 4, Ry, 4, 4, tmp );
  88.    copy_matrix( tmp, Tw, 4, 4 );
  89.  
  90.    set_to_identity( Td, 4 );
  91.    Td[3][1] = SCREEN_HEIGHT;
  92.    matrix_mult( Tw, 4, 4, Td, 4, 4, tmp );
  93.    copy_matrix( tmp, Tw, 4, 4 );
  94.    /* Now, Tw has the world to device/screen transformation matrix. */
  95.    if ( proj == PARALLEL )
  96.    {
  97.       /* Beautiful!!!  Perform a single transformation of the image (for
  98.          parallel projection). */
  99.       printf( "Transforming, projecting and mapping the image onto screen." );
  100.       matrix_mult( Tm, 4, 4, Tw, 4, 4, tmp );
  101.       copy_matrix( tmp, Tm, 4, 4 );
  102.       show_matrix( Tm, 4, 4 );
  103.       transform_and_homogenize_image( image, tr_image, Tm );
  104.       printf( "Done.\n" );
  105.       printf( "Rendering the image.  " );
  106.       render_image( LEFT_EYE );        printf( "Done.\n" );
  107.    }
  108.    if ( proj == PERSPECTIVE )
  109.    {
  110.       printf( "Transforming, projecting, homogenizing and mapping the " );
  111.       printf( "Left image onto screen.  " );
  112.       matrix_mult( Tm, 4, 4, Left, 4, 4, tmp );
  113.       matrix_mult( tmp, 4, 4, Tw, 4, 4, Left );
  114.       printf( "Left eye transformation matrix.\n" );
  115.       show_matrix( Left, 4, 4 );
  116.       transform_and_homogenize_image( image, tr_image, Left );
  117.       printf( "Done.\n" );
  118.       printf( "Rendering the Left eye image.  " );
  119.       render_image( LEFT_EYE );        printf( "Done.\n" );
  120.  
  121.       printf( "Transforming, projecting, homogenizing and mapping the " );
  122.       printf( "Right image onto screen.  " );
  123.       matrix_mult( Tm, 4, 4, Right, 4, 4, tmp );
  124.       matrix_mult( tmp, 4, 4, Tw, 4, 4, Right );
  125.       /* Move the right eye view into the lower half of the buffer. */
  126.       set_to_identity( Tw, 4 );
  127.       Tw[3][1] = SCREEN_HEIGHT;                 /* move in device coord */
  128.       matrix_mult( Right, 4, 4, Tw, 4, 4, tmp );
  129.       copy_matrix( tmp, Right, 4, 4 );
  130.       printf( "Right eye transformation matrix.\n" );
  131.       show_matrix( Right, 4, 4 );
  132.       transform_and_homogenize_image( image, tr_image, Right );
  133.       printf( "Done.\n" );
  134.       dszom( &0, &512, &1 );                   /* look at the lower half */
  135.       printf( "Rendering the Right eye image.  " );
  136.       render_image( RIGHT_EYE );        printf( "Done.\n" );
  137.    }
  138. }
  139.  
  140. [LISTING TWO]
  141.  
  142. make_composite_viewing_matrix( proj )
  143. int proj;      /* PARALLEL or PERSPECTIVE */
  144. {
  145.    double p1[4], p2[4], p3[4], p_tmp[4];
  146.    double T[4][4], Rx[4][4], Ry[4][4], Rz[4][4], C_tmp1[4][4], C_tmp2[4][4];
  147.    double Per[4][4], d1, d2, d12, cos_ang, sin_ang;
  148.    /* Initialize the three points */
  149.    p1[0] = eye_pt.x;  p1[1] = eye_pt.y;  p1[2] = eye_pt.z;  p1[3] = 1.0;
  150.    p2[0] = p2[1] = p2[2] = 0.0;  p2[3] = 1.0;
  151.    p3[0] = p1[0] + vup.x;  p3[1] = p1[1] + vup.y;  p3[2] = p1[2] + vup.z;
  152.                            p3[3] = 1.0;
  153.    /* Magnitude of vector p1->p2 */
  154.    d12 = sqrt( p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2] );
  155.    /* Create the translation matrix. */
  156.    set_to_identity( T, 4 );
  157.    T[3][0] = -p1[0];   T[3][1] = -p1[1];   T[3][2] = -p1[2];
  158.    /* Translate the three points p1, p2, and p3 to the origin. */
  159.    matrix_mult( p1, 1, 4, T, 4, 4, p_tmp );
  160.    copy_matrix( p_tmp, p1, 1, 4 );
  161.    matrix_mult( p2, 1, 4, T, 4, 4, p_tmp );
  162.    copy_matrix( p_tmp, p2, 1, 4 );
  163.    matrix_mult( p3, 1, 4, T, 4, 4, p_tmp );
  164.    copy_matrix( p_tmp, p3, 1, 4 );
  165.  
  166.    d1 = sqrt( p2[0] * p2[0] + p2[2] * p2[2] );      /* length of projection */
  167.    cos_ang = -p2[2] / d1;
  168.    sin_ang =  p2[0] / d1;
  169.    /* Create the rotation about Y-axis matrix. */
  170.    set_to_identity( Ry, 4 );
  171.    Ry[0][0] = cos_ang;  Ry[0][2] = -sin_ang;
  172.    Ry[2][0] = sin_ang;  Ry[2][2] =  cos_ang;
  173.    /* Rotate the three points p2, and p3 about the Y-axis. */
  174.    /* p1 is at the origin after translation => no need to rotate. */
  175.    matrix_mult( p2, 1, 4, Ry, 4, 4, p_tmp );
  176.    copy_matrix( p_tmp, p2, 1, 4 );
  177.    matrix_mult( p3, 1, 4, Ry, 4, 4, p_tmp );
  178.    copy_matrix( p_tmp, p3, 1, 4 );
  179.  
  180.    cos_ang = -p2[2] / d12;
  181.    sin_ang = -p2[1] / d12;
  182.  
  183.    /* Create the rotation about X-axis matrix. */
  184.    set_to_identity( Rx, 4 );
  185.    Rx[1][1] = cos_ang;   Rx[1][2] = sin_ang;
  186.    Rx[2][1] = -sin_ang;  Rx[2][2] = cos_ang;
  187.    /* Rotate the three points p2, and p3 about the X-axis. */
  188.    matrix_mult( p2, 1, 4, Rx, 4, 4, p_tmp );
  189.    copy_matrix( p_tmp, p2, 1, 4 );
  190.    matrix_mult( p3, 1, 4, Rx, 4, 4, p_tmp );
  191.    copy_matrix( p_tmp, p3, 1, 4 );
  192.    /* Sanity check. */
  193.    printf( "The view vector should be [ %lf %lf %lf %lf ]\n", 0.0, 0.0, -d12,
  194.                                                               1.0 );
  195.    printf( "The view vector is [ %lf %lf %lf %lf ]\n", p2[0], p2[1],
  196.                                                        p2[2], p2[3] );
  197.    d2 = sqrt( p3[0] * p3[0] + p3[1] * p3[1] );
  198.    cos_ang = p3[1] / d2;
  199.    sin_ang = p3[0] / d2;
  200.    /* Create the rotation about Z-axis matrix. */
  201.    set_to_identity( Rz, 4 );
  202.    Rz[0][0] = cos_ang;   Rz[0][1] = sin_ang;
  203.    Rz[1][0] = -sin_ang;  Rz[1][1] = cos_ang;
  204.    /* At this point the translation, and all rotation matrices are known
  205.       and need to be combined into a single transformaation matrix. */
  206.    matrix_mult( T, 4, 4, Ry, 4, 4, C_tmp1 );
  207.    matrix_mult( C_tmp1, 4, 4, Rx, 4, 4, C_tmp2 );
  208.    matrix_mult( C_tmp2, 4, 4, Rz, 4, 4, C_tmp1 );
  209.    copy_matrix( C_tmp1, Tm, 4, 4 );
  210. }
  211.  
  212.  
  213. [LISTING THREE]
  214.  
  215. void compute_normals()
  216. {
  217.    register i, j;
  218.    point_3D_t p11_p21, p11_p22, p11_p12, cross_a, cross_b;
  219.    double d, dx;              /* stepping distance (dx = dy) */
  220.    double dx_sqrd;            /* dx^2 */
  221.    int num_lines;
  222.    point_3D_t na, nb;         /* normals to triangle A and B */
  223.    dx = scan_sz / num_samples;
  224.  
  225.    num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  226.    num_lines--;
  227.    dx_sqrd = dx * dx;
  228.    for( i = 0; i < num_lines; i++ )
  229.       for( j = 0; j < num_lines; j++ )
  230.       {
  231. #if DEBUG
  232.          printf( "P11 %f %f %f\n", image[i][j].x, image[i][j].y,
  233.                                    image[i][j].z );
  234.          printf( "P21 %f %f %f\n", image[i+1][j].x, image[i+1][j].y,
  235.                                    image[i+1][j].z );
  236.          printf( "P12 %f %f %f\n", image[i][j+1].x, image[i][j+1].y,
  237.                                    image[i][j+1].z );
  238.          printf( "P22 %f %f %f\n", image[i+1][j+1].x, image[i+1][j+1].y,
  239.                                    image[i+1][j+1].z );
  240. #endif
  241.          p11_p21.z = image[ i + 1 ][ j ].z - image[ i ][ j ].z;
  242.          p11_p22.z = image[ i + 1 ][ j + 1 ].z - image[ i ][ j ].z;
  243.          p11_p12.z = image[ i ][ j + 1 ].z - image[ i ][ j ].z;
  244. #if DEBUG
  245.          printf( "dz11_21 = %f,  dz11_22 = %f,  dz11_12 = %f\n",
  246.                  p11_p21.z, p11_p22.z, p11_p12.z );
  247. #endif
  248.          /* It's possible to eliminate one more multiplication in the
  249.             computations below. */
  250.          na.x = dx * ( p11_p21.z - p11_p22.z );
  251.          na.y = dx * p11_p21.z;
  252.          na.z = dx_sqrd;
  253. #if DEBUG
  254.          printf( "Na %f %f %f\n", na.x, na.y, na.z );
  255. #endif
  256.          nb.x = (-dx) * p11_p12.z;
  257.          nb.y = dx * ( p11_p22.z - p11_p12.z );
  258.          nb.z = dx_sqrd;
  259. #if DEBUG
  260.          printf( "Nb %f %f %f\n", nb.x, nb.y, nb.z );
  261. #endif
  262.          /* Normalize the normal vectors, since the intensity will be
  263.             proportional to the angle between light source and the normal. */
  264.          d = sqrt((double)( na.x * na.x + na.y * na.y + na.z * na.z ));
  265.          na.x = na.x / d;
  266.          na.y = na.y / d;
  267.          na.z = na.z / d;
  268. #if DEBUG
  269.          printf( "Na %f %f %f\n", na.x, na.y, na.z );
  270. #endif
  271.          d = sqrt((double)( nb.x * nb.x + nb.y * nb.y + nb.z * nb.z ));
  272.          nb.x = nb.x / d;
  273.          nb.y = nb.y / d;
  274.          nb.z = nb.z / d;
  275. #if DEBUG
  276.          printf( "Nb %f %f %f\n", nb.x, nb.y, nb.z );
  277. #endif
  278.          /* Compute the dot product between the light source vector and
  279.             the normals (== to the angle between two unit vectors ).
  280.             -1 <= cos( theta ) <= 1, which will be very useful. */
  281.          image[ i ][ j ].sha = light_source.x * na.x + light_source.y * na.y +
  282.                                light_source.z * na.z;
  283.          image[ i ][ j ].shb = light_source.x * nb.x + light_source.y * nb.y +
  284.                                light_source.z * nb.z;
  285.       }
  286. }
  287.  
  288.  
  289. [LISTING FOUR]
  290.  
  291. transform_and_homogenize_image( s, d, tm )
  292. point_3D_ex_t  s[][ MAX_IMAGE_SIZE ], d[][ MAX_IMAGE_SIZE ];
  293. double *tm;         /* transformation matrix */
  294. {
  295.    register i, j;
  296.    int  num_lines;
  297.    double  p[4];      /* the point to be transformed */
  298.    double  t[4], inv_W;
  299.  
  300.    num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  301.    for( i = 0; i < num_lines; i++ )
  302.       for( j = 0; j < num_lines; j++ )
  303.       {
  304.          p[0] = s[i][j].x;   p[1] = s[i][j].y;
  305.          p[2] = s[i][j].z;   p[3] = 1.0;
  306.          matrix_mult( p, 1, 4, tm, 4, 4, t );
  307.          if ( t[3] != 1.0 )        /* divide by W (homogenize) */
  308.          {
  309.             inv_W = 1.0 / t[3];
  310.             t[0] *= inv_W;   t[1] *= inv_W;  t[2] *= inv_W;
  311.          }
  312.          d[i][j].x = t[0];   d[i][j].y = t[1];   d[i][j].z = t[2];
  313.          d[i][j].sha = s[i][j].sha;    d[i][j].shb = s[i][j].shb;
  314.       }
  315. }
  316.  
  317.  
  318. [LISTING FIVE]
  319.  
  320. render_image( y )
  321. int  y;
  322. {
  323.    register i, j;
  324.    int    num_lines;
  325.    short  v[6], intensity;
  326.  
  327.    num_lines = MIN( num_samples, MAX_IMAGE_SIZE );
  328.    num_lines--;
  329.    for( i = 0; i < num_lines; i++ )
  330.       for( j = 0; j < num_lines; j++ )
  331.       {
  332.          v[0] = ROUND( tr_image[ i     ][ j     ].x );
  333.          v[1] = ROUND( tr_image[ i     ][ j     ].y );
  334.          v[2] = ROUND( tr_image[ i + 1 ][ j     ].x );
  335.          v[3] = ROUND( tr_image[ i + 1 ][ j     ].y );
  336.          v[4] = ROUND( tr_image[ i + 1 ][ j + 1 ].x );
  337.          v[5] = ROUND( tr_image[ i + 1 ][ j + 1 ].y );
  338.          /* Render triangle A */
  339.          intensity = ROUND( tr_image[ i ][ j ].sha * (double)(NUM_SHADES - 1));
  340.          if ( intensity > ( NUM_SHADES - 1 ))
  341.             intensity = NUM_SHADES - 1;      /* saturate */
  342.          if ( intensity < 0 )
  343.          {
  344. #if 0
  345.             printf( "Triangle A,  intensity = %d\n", intensity );
  346.             printf( "v11.x = %f, v11.y = %f, v11.z = %f\n",
  347.                     image[i][j].x, image[i][j].y, image[i][j].z );
  348.             printf( "v21.x = %f, v21.y = %f, v21.z = %f\n",
  349.                     image[i+1][j].x, image[i+1][j].y, image[i+1][j].z );
  350.             printf( "v22.x = %f, v22.y = %f, v22.z = %f\n",
  351.                     image[i+1][j+1].x, image[i+1][j+1].y, image[i+1][j+1].z );
  352. #endif
  353.             intensity = 0;
  354.          }
  355.          if ( clip_to_viewport( v, 6, y ) == ACCEPT )
  356.             dspoly( &intensity, v, &6 );
  357.  
  358.          v[2] = ROUND( tr_image[ i     ][ j + 1 ].x );
  359.          v[3] = ROUND( tr_image[ i     ][ j + 1 ].y );
  360.          /* Render triangle B */
  361.          intensity = ROUND( tr_image[ i ][ j ].shb * (double)( NUM_SHADES-1));
  362.          if ( intensity > ( NUM_SHADES - 1 ))
  363.             intensity = NUM_SHADES - 1;      /* saturate */
  364.          if ( intensity < 0 )   intensity = 0;
  365.          if ( clip_to_viewport( v, 6, y ) == ACCEPT )
  366.             dspoly( &intensity, v, &6 );
  367.       }
  368. }
  369.  
  370.  
  371.  
  372. [LISTING SIX]
  373.  
  374. display_3D_data( proj )
  375. int proj;
  376. {
  377.    double  Tw[4][4], S[4][4], Td[4][4], Ry[4][4], Per[4][4], tmp[4][4];
  378.    double  Left[4][4], Right[4][4];
  379.    double  e_v,      /* interocular distance mapped back to the view  coord. */
  380.            e_w;      /* interocular distance mapped back to the world coord. */
  381.  
  382.    printf( "Computing normals for shading.  " );
  383.    compute_normals();      /* and the dot products for each triangle */
  384.    printf( "Done.\n" );
  385.  
  386.  
  387.  
  388.